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I  ABSTRACT 

An  efficient  method  is  described  for  sensitivity  analysis  of  nonlinear 
initial  value  problems,  which  may  Include  algebraic  equations  as  well  as 
ordinary  differential  equations .{ODE's ).CJ 

The  linearity  of  the  sensitivity  equations  is  utilized  to  solve  them 
directly  via  the  local  Jacobian  of  the  state  equations.  The  method  is 
Implemented  with  the  implicit  Integrator  DASSL  and  is  demonstrated  on  a 
stiff  Industrial  reaction  model. 


AMS  (MOS)  Subject  Classifications:  34-04,  65L05,  65H10 

Key  Words:  Numerical  integration,  Stiff  equations.  Parametric  sensitivities. 
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SIGNIFICANCE  AND  EXPLANATION 


Many  physical  systems  are  modelled  by  systems  of  ordinary  differential 
and  algebraic  equations  with  initial  conditions.  The  solution  vector  u 
depends  on  the  time  t,  and  on  a  vector  J3  of  unknown  parameters.  This 
report  deals  with  the  calculation  of  the  first-order  parametric  sensitivities, 
Wik(t,e)  =  Su./ae^,  which  are  useful  in  parameter  estimation,  system  design 
and  control . 

The  method  giv--n  here  takes  advantage  of  the  similarity  of  the  backward 
difference  forms  of  the  u-equations  and  W-equations,  as  well  as  the  linearity 
of  the  W-equatlons,  to  achieve  unusually  fast  solutions  with  minimal  memory 
requirements.  The  method  has  been  implemented  as  a  modification  of  the 
program  package  DASSL.  Numerical  results  are  given  for  a  chemical  kinetic 
example. 


The  responsibility  for  the  wording  and  views  expressed  in  this  descriptive 
summary  lies  with  MRC,  and  not  with  the  authors  cf  this  report. 


With  the  rapid  development  of  digital  computers.  Increasingly  realistic 
mathematical  models  are  being  used  to  investigate  chemical  phenomena.  New 
mechanistic  features,  however,  call  for  new  physicochemical  parameters  whose 
values  may  not  be  accurately  known.  Consequently,  there  Is  an  Increasing 
need  for  parametric  sensitivity  analysis  of  proposed  differential  and 
algebraic  models. 

Parametric  sensitivity  analysis  Is  a  very  active  research  area. 
Extensive  reviews  can  be  found  in  Rabitz  et  al.  [14],  and  In 
Tilden  et  al.  [17],  Applications  occur  In  every  engineering  and  scientific 
discipline.  Potential  areas  of  application  In  chemical  engineering  Include 
optimization,  parameter  estimation,  model  simplification,  process 
sensitivity  and  multiplicity,  experimental  design  and  many  more. 

In  this  paper  we  address  the  problem  of  numerical  computation  of 
sensitivity  functions  for  systems  of  ordinary  differential  and  algebraic 
equations.  We  develop  a  simple,  efficient  algorithm  for  this  purpose  by 
extension  of  a  standard  implicit  Integrator. 


jENoITU'ITY  ANALYSIS  OF  INITIAL  VALUE  PROBLEMS 
WITH  MIXED  ODE'S  AND  ALGEBRAIC  EQUATIONS 

Makis  Caracotsios  and  Warren  E.  Stewart 


PROBLEM  STATEMENT. 

Consider  a  dynamic  system,  described  by  the  following  set  of 
differential  and  algebraic  equations: 

EH'(t)  -  f(t,y(t);g)  (la) 

-  H0(fi)  (lb) 

Here  u  Is  an  n-dlmensional  vector  of  state  variables,  8  Is  an  m-dlmenslonal 

*  m 

vector  of  time-independent  parameters  and  E  Is  an  (n,n)  matrix  of  constant 

coefficients.  Most  frequently  In  chemical  kinetics  calculations  the  matrix 

E  assumes  the  form 
* 


E 


L2 


a 


(2) 


where  l's'  Is  the  Identity  matrix  of  order  s.  If  s  ■  n,  system  (la)  con¬ 
sists  of  purely  differential  equations.  If  1  <  s  <  n,  system  (la)  consists 
of  ordinary  differential  and  algebraic  equations.  The  latter  case  arises, 
for  Instance,  In  analysis  of  reaction  schemes  where  equilibria  give 
algebraic  constraints  on  the  concentrations. 

We  define  the  (n,m)  matrix  W(t)  of  sensitivity  functions  as 


¥(t): 


3«(t) 


(3) 
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This  matrix  satisfies  a  set  of  dlfferentlal/algebralc  equations  which  can 

be  derived  by  partial  differentiation  of  equations  (la),  (lb)  with  respect 
to  the  parameter  vector  6: 


En'(t)  -  J(t)t|(t)  ■  t(t,||(t>;8> 

38 


H<ft0) 


aU0(fl> 


(4a) 

(4b) 


where  the  matrix  jJ(t)  (shorthand  for  ij(t,u(t);§))  Is  defined  as 

ij(t):  ■  |j  f(t,y(t);g)  (5) 

Various  properties  of  equations  (4a)  and  (4b)  are  described  In  Tomovlc 
and  Vukobratovlc  [18].  The  most  striking  feature  of  these  sensitivity 
equations  Is  that  they  are  linear,  regardless  of  the  linearity  or 
nonlinearity  of  the  state  equations  (la)  and  (lb).  The  problem  studied  here 
Is  the  numerical  computation  of  the  matrix  W(t)  from  equations  (4a) 
and  (4b). 


LITERATURE  REVIEW  AND  THEORETICAL  BACKGROUND 

Before  describing  the  new  algorithm,  we  review  a  few  known  facts  about 
the  solution  of  mixed  systems  of  ordinary  differential  and  algebraic 
equations.  Several  Investigators  [6],  [73,  [15],  have  considered  this 
subject  and  recently  Petzold  [11]  has  published  an  algorithm  called  DASSL 
for  the  solution  of  such  systems. 

Not  all  systems  of  dlfferentlal/algebralc  equations  are  solvable.  The 
reader  Is  encouraged  to  consult  the  literature,  Petzold  and  Gear  [12],  and 
Campbell  and  Petzold  [2],  on  this  peculiar  feature  of  mixed  systems. 
However,  for  the  systems  that  we  are  considering,  where  the  matrix 
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E  assumes  the  form  (2),  sufficient  conditions  are  known  [12]  for  the 
solvability  of  (la)  and  (lb). 


Let  the  function  f(t,u(t);0)  be  continuously  differentiable 
with  respect  to  u(t).  Now  consider  the  Jacobian  matrix  J(t)  defined  by 
equation  (5)  for  the  system  (la).  If  we  partition  the  Jacobian  matr'/ 
according  to  the  partition  of  E,  l.e., 


J(t) 


Jl2(t»‘ 

^21^ 

^22^ 

(6) 


where  Ju(t)  Is  an  (s,s)  matrix,  then  the  system  (la)  and  (lb)  Is 
solvable  If 


det  jj22(t)  *  0  for  a11  *  (7) 

Under  this  condition,  the  solution  obtained  by  a  k-step  backward 
differentiation  formula  algorithm  with  k  '  7  and  fixed  step  size  h  converges 
to  0(hk)  If  all  Initial  values  are  correct.  Further  aspects  of 
equations  (la)  and  (lb)  and  their  numerical  treatment  are  discussed 
In  Petzold  [13]. 

Let  us  now  review  some  of  the  methods  used  for  the  computation  of  the 
sensitivity  matrix  W(t).  With  one  exception  (Stewart  and  Sorensen  [16]) 
the  known  methods  are  for  systems  of  ODE 1 s  only;  that  Is,  for  systems  with 
E  =  I'"'.  The  available  algorithms  Include  the  Fourier  amplitude  test  [4], 
direct  differential  methods  [5],  Green's  function  methods  [8],  the 
analytically  Integrated  Magnus  method  [14]  (a  modification  of  the  Green's 
function  method),  and  finite  difference  methods.  The  Green's  function 


method  and  Its  variations  exploit  the  fact  that  the  sensitivity  equations 
are  linear  inhomogeneous  with  time  varying  coefficients;  consequently  they 
can  be  solved  by  first  calculating  the  solution  of  the  homogeneous  part  and 
then  determining  the  particular  solution  corresponding  to  each  parameter. 

Several  authors  have  proposed  to  solve  the  sensitivity  equations  by 
extending  known  solution  algorithms  for  the  state  equations.  This  idea  is 
based  on  the  Identity  of  the  coefficient  matrices  in  the  sensitivity 
equations  to  those  in  the  locally  linearized  form  of  the  state  equations 
on  the  locus  u(t).  Stewart  and  Sorensen  [16],  Vemuri  and  Raefsky  [19], 

Lojek  [10],  and  Hwang  [9]  have  developed  various  aspects  of  this  method; 
nevertheless  the  idea  is  still  under  development. 

In  the  present  work,  we  exploit  fully  the  similarity  of  the  sensitivity 
and  state  equations,  by  building  the  sensitivity  analysis  into  a  robust 
differential/algebraic  equation  solver.  Then  we  Illustrate  the  algorithm 
by  solving  a  stiff  industrial  kinetics  problem. 


MATHEMATICAL  DEVELOPMENT  OF  THE  SENSITIVITY  ANALYSIS  ALGORITHM 


One  of  the  most  Important  steps  In  developing  the  sensitivity  analysis 
algorithm  was  the  selection  of  the  Integrator.  For  mixed  systems  of 
differential  and  algebraic  equations,  there  are  several  codes 
[7],  [11],  [15]  designed  to  perform  the  integration.  These  codes  are 
primarily  based  on  an  idea  developed  by  Gear  [6];  specifically,  the 
derivative  u'(tn+1)  is  approximated  by  a  backward  difference  formula  with 
adaptable  order  and  step  size,  and  the  resulting  system  of  nonlinear 
equations  is  solved  for  u( tn+ ^ )  via  a  modified  Newton  scheme. 

We  chose  for  this  work  the  package  DASSL  developed  by  Petzold  [11]. 


This  package  Is  portable,  robust  and  easy  to  use.  We  tested  the  code 
successfully  on  a  wide  variety  of  stiff  problems,  both  differential  and 
mixed,  before  adding  the  sensitivity  analysis  algorithm. 

First,  consider  the  solution  of  state  and  sensitivity  equations 
(la),  (lb),  (4a),  (4b)  as  a  single  system.  In  this  approach,  one  needs  the 
Jacobian  matrix  of  the  total  system  (la)  and  (4a).  If  we  partition  the 
sensitivity  matrix  W(t)  Into  column  vectors  as 

y(t)  -  [^(t)  lw2(t)l 

where 

3y(t)  ... 

#,(t):  -  — —  1  B  1,2,  •••  m,  (9) 

*  o6j 

then  the  Jacobian  matrix  J*(t)  (shorthand  for  J*(t,u( t) »W( t) ;9) ) 
of  the  total  system  (la}  and  (4a}  Is 


git)  a  a . a 

jjiti  git)  a . a 

j2(t)  a  g<t>  a 

g*(t)  ■  .  ,  . 

•  •  •  • 

•  •  •  • 

•  •  •  • 

4m<t)  a  a . git) 


where 

3J( t)  3^(t) 


(10) 


(11) 
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The  evaluation  of  J*(t)  Is  a  formidable  calculation,  though  a  natural 
requirement  of  Newton  Iteration  on  the  total  equation  system,  (la)  and  (4a). 
A  simpler  and  quicker  approach  Is  to  solve  (la)  before  (4a)  at  each  time 
step,  as  shown  below;  then  the  matrices  j^(t)  are  not  required. 

Let  u^(t)  be  the  local  Interpolant  of  u(t)  obtained  in  r  Newton 
Iterations  of  a  kth-order  Integrator  within  a  given  time  step.  Then  the 
next  Iteration  will  give  the  Interpolant  u^(t)  +  Au^(t),  which 
satisfies  the  following  linearized  form  of  equation  (la) 


k(j5  (r*(t)+£y  (r)(t))  *  J(t,^r^(t);§)+jJ(t,j5^r^(t);§)^y^r)  +  0(hk) 
Hence  the  correction  satisfies 


(12) 


Wa'(r)N(r)(t)^5(r)  ■  t,y (r) <t) ;Q)  -  gj'(r)(t)  +  0(hk)  (13) 

when  the  standard  Newton  method  (with  J^(t)  updated  for  each  Iteration) 
Is  used.  If  Au^r)(t)  converges  toward  zero  with  Increasing  r,  then 
J^(t)  converges  toward  J(t),  and  equation  (13)  becomes  formally  similar 
to  equation  (4a).  Therefore,  we  can  defer  consideration  of  equation  (4a) 

f*  /  M  \  •» 

until  u'  (t)  has  converged  to  u(t)  at  the  current  value  of  t.  Then  we 
can  update  the  sensitivity  solution  W( t )  directly  by  use  of  equation  (4a), 
which  has  the  same  coefficients  as  equation  (13)  but  a  different,  now 
computable  right  hand  function.  More  specifically  the  corrections 
AWj(t)  are  calculated  via  a  single  iteration  by  solving 

Mii<tH(t)  ^(t)  =  <p)  (t)+J(  t)  i|p)(t)  +  J(t,y(t);§)  +  0(hk) 

1  =  1,2,  •••  m  (14) 
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in  which  wjP^(t)  is  the  predicted  value  of  w^(t)  via  a  kth-order 
predictor  formula. 

The  vectors  5j.f(t)  are  calculated  at  the  current  t  as  follows: 

tj.|(t)  =  +  ^(t)  i  =  1,2,  •••  m  (15) 

On  completion  of  the  update,  the  local  truncation  error  is  tested,  and  if 
necessary  the  step  size  h  and  approximation  order  k  are  adjusted  to  achieve 
the  specified  accuracy  for  u(t)  and  W(t). 

Numerical  tests  show  that  stringent  tolerances  on  u(t)  normally  lead 
to  a  good  solution  for  W(t)  as  well,  provided  that  J(t)  Is  updated  before 
computing  AW ( t ) .  On  the  other  hand,  If  the  Iteration  matrix  Is  only 
occasionally  updated  (as  Is  usual  In  Implicit  Integrators),  then  the  local 
tolerances  on  W(t)  are  essential  to  control  the  calculation. 

COMPUTER  IMPLEMENTATION 

The  computer  Implementation  of  the  sensitivity  analysis  algorithm  was 
done  as  follows: 

1)  The  working  arrays  used  by  DASSL  were  modified  to  provide  storage 
allocation  for  the  sensitivity  functions. 

2)  An  algorithm  was  written  for  the  automatic  formation  of  the 
sensitivity  equations. 

3)  The  Integrator  In  DASSL  was  properly  extended  to  Include  the 
solution  of  the  sensitivities. 

The  following  definitions  were  adopted: 


1)  Iteration  matrix  D(t):  =  cE  -  J(t),  where  c  is  a  constant  which 
depends  on  step  size  history 

2)  Residual  of  the  state  equations 

g0(t) :  =  (t)  -  f(t,y(t);g) 

3)  Residual  of  the  sensitivity  equations 

(t)  -  jJ(t)^{t)  -  |(t,y(t);§)  1  =  1,2,  •••  m 

The  calculation  sequence  is  outlined  in  the  adjoining  flow  chart. 
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Update  u(t),  RQ(t) 


Is  it^ntion  matrix  currant 


NUMERICAL  EXAMPLE 
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i.\ 


v- 
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The  algorithm  was  tested  on  a  wide  variety  of  problems,  ranging  from 
linear  and  nonlinear  systems  of  differential  equations  to  large  systems  of 
differential  and  algebraic  equations.  We  present  here  a  batch  reactor 
example  given  by  the  Dow  Chemical  Company  [1].  Figure  1  shows  the 
proposed  mechanism  for  the  reaction  system.  The  time  dependent 
concentrations  are  modelled  by  the  following  system  of  differential  and 
algebraic  equations: 


u ^ ( t )  -  ~k2u2( t)Ug(t) 


U2(t) 


u3(t) 


*  -k1u2(t)u6(t)  +  k_1u10( t)  -  k2u2(t)u8(t) 
■  k2u2(t)u8(t)  +  k3u4( t)u6(t)  -  k_3ug(t) 


u4(t) 


Ug(t) 


u6(l) 

U7(t) 


U8(t) 


u9<t) 


*  -k3u4(t)u6(t)  +  k_3Ug(t) 

-  k1u2(-t)u6(t)  -  k_iu10(t) 

■  -kjU^tlUgft)  -  k3u4(t)Ug(t)  +  k_jU10(t)  +  k_3u9(t) 

■  -0.0131  +  u6(t)  +  Ug(t)  +  u  g( t )  +  u^Ct) 

K2Ui(t) 

K2+u7(t) 

K3u3(t) 


(16) 


K3+u7(t) 


u 


K  iu  5  ( t ) 
kj+u  7'(tj 
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Figure  1.  Chemical  reaction  model  for  the  numerical  example  [1].  The 
numhers  are  used  for  the  enumeration  of  the  chemical  species 
as  in  equation  [16). 
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with  initial  conditions 


u2 ( 0)  =  1.5776 
u2(0)  =  8.32 

Uj(0)  =0  j  =  3,4,5,9,10 
ug { 0)  =  0.0131 

u7(0)  -  0.5{-K2  +/ K^MKgU^O)  } 
u8(0)  -  u y ( 0 ) 

The  following  values  of  rate  and  equilibrium  constants  were  used  [3] 


■  21.893 

hr”1  Kg  gmole"1 

-  2.14  E09 

hr”1 

*  32.318 

hr”1  Kg  gmole”1 

CO 

CO 

CVJ 

a 

hr”1  Kg  gmole”1 

«*  1.07  E09 

hr”1 

■  7.65  E— 18 

gmole  Kg’1 

=  4.03  E-ll 

gmole  Kg"1 

=  5.32  E-18 

gmole  Kg"1 

The  natural  logarithms  of  these  constants  make  up  the  parameter  vector  0. 

3u(t) 

Our  goal  Is  to  estimate  the  sensitivity  matrix  W(t)  =  — - —  . 

-  36 

The  combined  system  of  state  and  sensitivity  functions  consists  of 
90  equations,  54  of  which  are  differential  and  36  are  algebraic.  This 
problem  presents  a  severe  test  for  the  DASSL  integrator  and  the  sensitivity 
analysis  algorithm. 

Table  1.  summarizes  the  computational  effort  for  the  solution  of  the 
above  problem  on  a  VAX  11/780  computer.  All  calculations  were  carried  out 
In  double  precision,  A  mixed  local  truncation  error  control  provided  in 
[11]  was  used.  The  tolerances  for  the  sensitivities  were  equated  to  the 
tolerances  of  the  state  variables.  The  total  reaction  time  considered 
was  53  hr. 

Figures  2  and  3  show  the  evolution  of  the  concentration  profiles  as 
a  function  of  time,  while  Figures  4  through  11  show  the  dynamic  behavior 
of  the  sensitivity  functions. 

Sensitivity  plots,  like  those  In  Figures  4-11,  can  be  of  considerable 
use  to  the  theoretician  as  well  as  to  the  experimentalist.  From  Figures 
4  through  11  we  see  that  all  of  the  rate  and  equilibrium  constants  have 
comparable  effects  on  the  concentrations  and  therefore  we  cannot  eliminate 
any  step  from  the  proposed  mechanism  In  Figure  1.  However,  a  close 
Inspection  of  the  results  reveals  the  following  linear  relations: 


3y(t)  3u(t) 

k  .  -  -  -  K,  J - 

_1  sk^  l  3K1 

3y(t)  3u(t) 

k  i  — -  =  K„ - 


(19) 


(20) 


TABLE  1 


Computational  effort  for  the  solution  of  the  Dow  problem 

STATE  EQUATIONS  SENSITIVITY  EQUATIONS 


e*-1.0E-6 

e-l.OE-7 

e-l.OE-6 

e-l.OE-7 

Time  steps  : 

191 

265 

187 

264 

Function  evaluations: 

418 

587 

383 

534 

CPUsecs  : 

4.0 

5.5 

23.5 

33 

local  error  test  at  each  time  step  which  requires  roughly  that 
abs  (local  error)  less  or  equal  than  rtol »abs(u)+atol 
where  rtol =a to! “e. 


Figure  2 


Figure  2.  Concentration  functions  u^(t)  in  gmole/kg 
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Figure  7 


Figure  7.  Elements  of  the  sensitivity  vector  3u(t)/3£nk, 
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Figure  11. 

.  Elements  of  the  sensitivity  vector  3u(t)/aa.nk 
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Thus,  one  cannot  estimate  all  parameters  of  the  model  from  data  on  u(t). 
Equations  (19)  and  (20)  Indicate  that  (k.^/Kj)  and  k_3K3  should 
be  used  as  parameters;  then  the  predicted  u(t)  Is  Invariant  to  Kj  and  «3. 
This  Information  proved  useful  In  fitting  the  model  (16),  (17)  to  the 
published  experimental  data  for  this  system  [3], 

CONCLUSIONS  AND  SIGNIFICANCE 

An  efficient  method  has  been  developed  here  for  the  calculation  of 
parametric  sensitivities  for  mixed  systems  of  ordinary  differential  and 
algebraic  equations.  The  method  has  been  tested  successfully  on  many 
problems,  ranging  from  simple  ODE's  to  large  stiff  systems  of  differential 
and  algebraic  equations.  The  scheme  has  been  Implemented  for  a  robust 
Implicit  Integrator  (DASSL),  and  Is  applicable  to  any  Implicit  Integrator. 
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